The present invention relates to a method of mapping proton transverse relaxation time constants, or functions thereof, in a target subject to localized movement, such as abdominal tissue, using nuclear magnetic resonance imaging.
Magnetic resonance imaging is an imaging modality that has been developed to elucidate the internal structure of essentially diamagnetic bodies through exploitation of the phenomenon of nuclear magnetic resonance. From a semi-classical standpoint, in a static magnetic field, nuclear isotopes that have a nuclear magnetic moment experience a torque which causes the moments to precess around the axis of the field at a frequency that is proportional to the magnitude of the magnetic moment and the magnitude of the applied field. Further, the orientations of the nuclear magnetic moments are quantised in a limited number of spin states. Thus, for a particular nuclear species, the moments or spins precess at the same frequency in random phase around the direction of the field but in different equilibrium populations of spin states. If a radio frequency (RF) pulse is applied with a frequency that matches the precessional frequency of a certain nuclear species, the populations of the spin states will be perturbed from their equilibrium values. Further, the spins acquire a certain level of phase coherence in that they precess in some measure of synchrony with each other. After the RF pulse is removed, the spins return to their equilibrium populations by two relaxation processes for which a magnetic resonance signal having the same frequency as the RF pulse can be detected. One of the relaxation processes involves the return of the spins to their equilibrium population values, called spin-lattice or longitudinal relaxation, for which the relaxation rate is characterised by the longitudinal relaxation time constant T1. The other relaxation process is one in which the spins lose their phase coherence, called spin-spin or transverse relaxation, for which the relaxation rate is characterised by the transverse relaxation time constant T2. For a particular nuclear species, the relaxation rates can vary greatly according to the chemical environment surrounding each isotope, on both the molecular and macro-molecular scale. The outstanding image contrast that can be achieved by MRI is a function of the variation in these relaxation rates, coupled with the variations in nuclear density that occur throughout the body being imaged.
Although there are a number of nuclear species for which magnetic resonance can be observed, the hydrogen proton is of the greatest relative sensitivity. Consequently, it is the nuclear species around which magnetic resonance imaging has been developed. The hydrogen proton is also the most abundant nuclear species within the human body, with approximately two thirds of the body hydrogen contained in water molecules and the remainder found in fat and protein. The hydrogen proton thus makes an ideal probe for anatomical imaging. The remarkable level of soft tissue contrast that can be obtained by MRI is a result of the variation in hydrogen proton density and relaxation times for different tissues, and the perturbation of these times in various disease states.
Magnetic resonance images are constructed by varying the magnetic field strength in three dimensions throughout the subject or target to be examined. The variations in magnetic field result in precessional frequency changes of the nuclear species at various points in space, and thus enable the discrimination of magnetic resonance signals from different spatial locations. A map of signal intensities can then be constructed to obtain a magnetic resonance image. Depending on the manner in which the RF pulse is applied to cause magnetic resonance, the images that result can either be predominantly T1-weighted or T2-weighted. In T1-weighted images, the image intensities predominantly reflect the progression of spin-lattice relaxation, the extent of which depends on when the magnetic resonance signals are acquired. In T2-weighted images, the intensities essentially reflect the progression of spin-spin relaxation.
Typically, the relaxation time constant for a given region of interest over the image sequence is determined through the fitting of an equation to the measured signal intensities that describes the return of the hydrogen protons to their equilibrium spin states. For relaxation processes modeled in this fashion, the relaxation time constant that is determined is essentially an average of each and every relaxation time constant for each and every hydrogen proton within the region of interest. However, within any given region of interest, there may be particular populations of hydrogen protons that do not necessarily neighbour in space but which neighbour in terms of the chemical and physical environments which the protons experience, and which are thus characterised by their own distinct relaxation times. For example, the population of hydrogen protons found in fat will have distinctly different relaxation times from the population of hydrogen protons found in extra-cellular water. Thus, depending on the number of images acquired at different measurement times, a number of relaxation processes may be determined within the one region of interest for different populations of hydrogen protons. For the number of relaxation processes that are desired to be resolved, the equation that describes the return of the hydrogen protons to their equilibrium spin states is separately summed for each distinct population. For transverse relaxation, where the average relaxation process is characterised by a single exponential decay term involving the relaxation time constant and the measurement time, the relaxation time constant calculated is typically referred to as that for single (or mono-) exponential decay. When two or more transverse relaxation processes are being modeled, a number of exponential decay terms are summed, and the resulting equation is referred to as one of multi-exponential decay. When only two transverse relaxation processes are being modeled, the equation is one for double (or bi-) exponential decay. In this instance, it is common to refer to fast and slow relaxation components of hydrogen protons, ie: a population of hydrogen protons that undergo fast relaxation back to their equilibrium spin rates, and a population of hydrogen protons that experience slow relaxation.
For a sequence of either T1- or T2-weighted images acquired at different measurement times, the relaxation time constants of the dominant relaxation process can be theoretically determined over the entire image. However, the calculation of a map of T1 or T2 relaxation times constants is complicated when the region to be examined is affected by some form of localized movement. This arises, for example, in imaging of the abdomen, where the regular, repetitive motion of breathing results in image intensity perturbations across the image. The existence of breathing artefacts over the region of interest makes the calculation of both accurate and complete T1 or T2 maps infeasible. To date, the successful generation of relaxation time maps has only been reported in those cases where sample movement is not a factor, as in non-medical applications of materials research and NMR microscopy, and in imaging of the brain, where pulsatile and respiratory affects may be considered to be negligible. The successful generation of accurate and complete relaxation time maps over the abdomen or other targets subject to a similar extent of localized movement has not been demonstrated.
The image intensity perturbations caused by localized movement arise as a result of the method of image construction. To obtain magnetic resonance images of sufficient intensity for both qualitative and quantitative analysis, the image signal intensities must be sufficiently above the background image noise. To obtain such intensities, repeated measurements must be performed over the same measurement (or repetition) time and the signal intensities cumulated. As noise is not additive, the signal-to-noise ratio increases, and useful image intensities can be obtained. However, this process is adversely affected when the area to be examined is subject to localized movement, as in imaging of the abdomen. The action of breathing causes the image signal intensity measured for a particular volume in space to be some average of the relaxation processes and proton densities for those parts of the body moving through that volume. Consequently, the image intensity no longer identifies a fixed volume within the subject. More importantly though, for a sequence of images acquired at different measurement times, the calculation of a map of accurate and complete relaxation time constants for some region of interest over the images becomes infeasible.
Compensating for movement within the magnetic resonance images does not in itself ensure the calculation of accurate relaxation time constants. To obtain truly valid relaxation time constants, a number of other factors must be accounted for.
It is an objective of the present invention to provide a method of generating a map or distribution of values of parameters that are a function of the proton transverse relaxation time constants in a target subject to localized movement, using nuclear magnetic resonance imaging.
According to the present invention there is provided a method of generating a map or distribution of values of parameters that are a function of proton transverse relaxation time constants T2 in a target subject to localized movement using nuclear magnetic resonance imaging, the method comprising at least the steps of:
acquiring a sequence of three or more T2-weighted magnetic resonance images of an identical image plane through the target at different spin-echo times;
defining a region of interest (ROI) over the target within the T2-weighted images that identifies the same region across all of the images;
processing image intensities of each T2 weighted image within the region of interest (ROI) by application of an image filter that compensates for image intensity perturbations caused by localized movement of the target in said ROI to produce a filtered image intensity for each T2-weighted image;
calculating parameters that are a function of T2 by a curve-fitting procedure applied to the filtered image intensities as a function of the spin-echo time (TE); and,
generating a map or distribution of, values of said parameters calculated by said curve fitting procedure or of, other parameters that are a function of or correlate with the parameters calculated by said curve fitting procedure.
Preferably said step of acquiring the sequence of T2-weighted images is performed at spin-echo times that are at least one order of magnitude less than the time between repetitions of consecutive RF pulse sequences used in the generation of said T2-weighted images. This time is known as the xe2x80x9crepetition timexe2x80x9d.
Preferably said application of an image filter involves employing a procedure of spatial neighbourhood averaging to replace an image intensity within any T2-weighted image at a given location within said ROI by an intensity that is a function of neighboring intensities in that image.
Preferably said neighboring intensities employed in the spatial neighbourhood averaging are bounded by a rectangular window kernel that covers a first range of movement of the target along a first axis of the image plane and a second range of movement of the target along a second perpendicular axis of the image plane.
Preferably, when said target is abdominal tissue, the first range of movement is between about 5 mm to 17 mm along said first axis perpendicular to a coronal plane and the second range of movement is about 4 mm to 14 mm along said second axis lies either in a sagittal plane or an axial plane whereby, in use, the method accommodates for the breathing artefacts in the T2-weighted images.
Preferably said curve fitting procedure includes incorporation of factors to compensate for one or more of background signal level offset; instrumental drift; and measurement errors on the image intensities. Preferably said curve fitting procedure also includes consideration of the vicinity of the image intensities to background noise level.
Preferably said curve fitting procedure includes fitting the following equation to the processed image intensities as a curve which models the decay of said intensities with increasing TE       I    ⁡          (      TE      )        =                    ∑                  n          =          1                N            ⁢              (                                            I              n                        ⁡                          (              0              )                                ⁢                      xe2x80x83                    ⁢                      ⅇ                                          -                TE                            /                              T                                  2                  ⁢                  n                                                                    )              +          S      LO      
where:
I(TE) is a processed image intensity at a given TE;
SLO is the signal level offset of the image signal intensities; and, when N=1, the said equation models single exponential decay, and where substituting I(0) for I1(0) and T2 for T21:
I(0) is the unknown intensity at TE=0 ms and is to be determined by the curve fitting procedure;
T2 is the unknown transverse relaxation time constant that characterises the decay of the processed image intensities with increasing TE and is to be determined by the curve fitting procedure;
or, when N=2, the said equation models bi-exponential decay, and where substituting If(0) for I1(0), T2f for T21, Is(0) for I2(0), and T2s for T22:
If(0) is the unknown intensity at TE=0 ms due to a fast relaxation component of hydrogen protons and is to be determined by the curve fitting procedure;
T2f is the unknown transverse relaxation time constant that characterises the decay of the processed image intensities with increasing TE for the fast relaxation component of hydrogen protons and is to be determined by the curve fitting procedure;
Is(0) is the unknown intensity at TE=0 ms due to a slow relaxation component of hydrogen protons and is to be determined by the curve fitting procedure;
T2s is the unknown transverse relaxation time constant that characterises the decay of the processed image intensities with increasing TE for the slow relaxation component of hydrogen protons and is to be determined by the curve fitting procedure.
Preferably said signal level offset is determined by analysis of background image intensities in a region free of image intensity perturbations caused by localized movement.
Preferably said background is a region within the T2-weighted images that is predominantly free of the presence of hydrogen protons.
Preferably said signal level offset is taken as the mean intensity of a Poisson distribution fitted to the distribution of background image intensities.
Preferably background noise level is derived as the mean of said Poisson distribution plus one standard deviation.
Preferably said Poisson distribution fitted to the background image intensities is the generalised Poisson distribution of the form             P      ⁡              (        χ        )              =                  α        ⁢                  xe2x80x83                ⁢                              θ            ⁡                          (                              θ                +                                  χ                  ⁢                                      xe2x80x83                                    ⁢                  λ                                            )                                            (                          χ              -              1                        )                          ⁢                  ⅇ                      -                          (                              θ                +                                  χ                  ⁢                                      xe2x80x83                                    ⁢                  λ                                            )                                                                        2            ⁢                          xe2x80x83                        ⁢            π            ⁢                          xe2x80x83                        ⁢            χ                          ⁢                              (                          χ              /              ⅇ                        )                    χ                      ,      xe2x80x83    ⁢      χ    ≥    1  
where:
"khgr" is an image intensity value;
xcex1 is a scale factor, and;
xcex8 and xcex parameterise the Poisson distribution,
and for which the mean xcexc of the Poisson distribution is
xcexc=xcex8(1xe2x88x92xcex)xe2x88x921
and the variance "sgr"2 is
"sgr"2=xcex8(1xe2x88x92xcex)xe2x88x923
Preferably the method further includes imaging of a phantom in-situ with the subject in the T2-weighted magnetic resonance images for image intensity correction or image intensity error calculation purposes.
Preferably the T2-weighted images are acquired with fixed gain settings.
Preferably said T2-weighted images are acquired over spin-echo times for which image intensities of the phantom are substantially constant, or for which a percentage change in the image intensities of the phantom over the spin-echo times is accurately known.
Preferably said phantom image intensities over the sequence of T2-weighted images are analysed to determine ROI scale factors for each image relative to one of the images to correct for the possible instrumental drift that perturb the image signal intensities in the ROI over the sequence of T2-weighted images.
Preferably said ROI scale factors that correct for the instrumental drift are determined by:
choosing one of the T2-weighted images as a reference image and treating any remaining T2-weighted images as scalable images;
obtaining mean intensity values for a number of regions over the phantom which are substantially free of intensity gradients in three or more of the T2-weighted images and for which the regions are the same in each image;
assigning the mean intensity values for the reference image to reference intensity values;
dividing the mean intensity value for every region over the phantom for each scalable image by the reference intensity value for a matching region in the reference image to obtain an instrumental drift scale factor for every region in each of the scalable images;
calculating the mean instrumental drift scale factor from the instrumental drift scale factors from any number of regions over the phantom for each scalable image;
using the standard deviation of the instrumental drift scale factors as an uncertainty factor on the mean instrumental drift scale factor, and;
assigning the mean instrumental drift scale factor for a scalable image to the ROI scale factor for that image.
Preferably the method also includes the step of incorporating any percentage change of said phantom image intensities over the spin-echo times into the ROI scale factor for each image for those T2-weighted image sequences for which the phantom image intensities are not effectively constant.
Preferably said instrumental drift over said sequence of T2-weighted images is corrected for by dividing the filtered image intensities of each image by an associated ROI scale factor for that image.
Preferably the regions over the phantom for which instrumental drift scale factors are calculated exhibit a substantially identical measure of RF field variation as that over the ROI.
Preferably said intensity measurement error on the filtered image intensities specified for the curve-fitting procedure is calculated as the standard error over the ROI of differenced intensity values which are the filtered image intensities subtracted from unprocessed image intensities.
Preferably said standard error used as the intensity measurement error on the filtered image intensities is calculated as the standard deviation over the ROI of said differenced intensity values divided by the square root of the number of intensity values in the neighbourhood of intensities over which the image filtering was performed.
Preferably said parameters calculated by said curve-fitting procedure for the ROI facilitate identification and characterisation of normal and abnormal tissue types.
Preferably said ROI is through the liver of an animal or person and said parameters calculated by the curve-fitting procedure are used to determine the presence and extent of one or more of: fibrosis, cirrhosis, lesions or tumours.
Preferably said parameters calculated by the curve-fitting procedure are formulated so that they further correlate with a measure of tissue iron overload for that region.
Preferably said formulation for when N=1 for the curve-fitting equation is:   1      T    2  
Preferably said formulation for when N=2 for the curve-fitting equation is:                     I        f            ⁡              (        0        )              -          S      LO                  (                                    I            f                    ⁡                      (            0            )                          +                              I            s                    ⁡                      (            0            )                          -                  2          ⁢                      S            LO                              )        ⁢          T              2        ⁢        f            
Preferably said measure of tissue iron overload is a measure of hepatic iron concentration.
Preferably said parameters calculated by said curve-fitting procedure are transverse relaxation rates and for which the distribution of transverse relaxation rates is parameterised by as few as one or two Gaussian functions to further characterise the ROI.